Agron508_assignment1

1.Equation to get PPFD

PPFD is photosynthetic photon flux density in mol/m2/s Irradiance is the radiant flux received by a surface per unit of area in W/m2 Energy content in the PAR waveband = 2.35*10^5 J/mol

# irr = irradiance
PPFD<-function(irr){   # here the function called PPFD used the irradiance "irr" as an input
    result<-irr*1/(2.35*10^5)  # result represents the output (that is the PPFD)
return(result) # you need to use return to be sure the result will be printed in your R console
}

PPFD(500)  # as PPFD function was created now you can put any input (irradiance)

## [1] 0.00212766

# if irr = 500 W/m2, PPFD is equal to 2100 micromol/m2/s

1.bis. PPFD graph

values<-seq(0,1000,by=100) 
# values is a vector including value of irradiance fromm 0 to 1000 with a increment of 100

ggplot(data=NULL,aes(x=values,y=PPFD(values)))+geom_point() +
  ylab("PPFD (mol/m2/s)") +
  xlab("Irradiance (W/m2)") +
  ggtitle("PPFD as a function of irradiance")

2.Planck’s equation function

Planck’s law describes the spectral density of electromagnetic radiation emitted by a black body in thermal equilibrium at a given temperature T

h<-6.626*10^-34 # Planck's constant (J*s)
c<-2.997925*10^8 # speed of light (m/s)  or 10^17 nm/s
k<-1.381*10^-23 # Boltzmann's constant (J/K)

planck<-function(wl,temp){  # wl = wavelength (m), T = temperature (K)
  num<-2*pi*h*c^2
  den<-wl^5*exp((h*c/k+wl+temp)-1)
  result<-num/den
  return(result)
}

3.Planck equation graph

h=6.62606896e-34
c=2.99792458e8
k=1.3806504e-23 
lambda = seq(1e-8,1e-4, by=1e-8)
T = 6000


y=(2*pi*h*(c^2))/(lambda^5)*(1./(exp((h*c)/(lambda*k*T))-1))


T2 = 288
y2=(2*pi*h*(c^2))/(lambda^5)*(1./(exp((h*c)/(lambda*k*T2))-1))
plot(1e6*lambda,1e-6*1e-6*y, xlim=c(0.0,5.0),ylim=c(0.0,100.0), xlab='Wavelength (um)', ylab='spectral emittance (MW/m2/um)')

# par(new=T)
# plot(1e6*lambda,1e-6*1e-6*y2, xlim=c(0.0,10.0), ylim=c(0.0,1/10000.0), xlab='', ylab='', axes=F)
# par(new=F)
plot(1e6*lambda,1e-6*1e-6*y2, xlim=c(0.0,50.0), ylim=c(0.0,3*1e-5), xlab='Wavelength (um)', ylab='spectral emittance (MW/m2/um)')

4.Active example Dark vision

4.1 Description step by step for one value of wavelength

lambda<-0.505 # micro m (or 505 nm)
# area pupil = 2.83*10^-5 m2
# 4.00*10^-11 W/m2 is the intensity for dark- adapted vision for human
# intensity per area = 4.00*10^-11 * 2.83*10^-5 = 1.13 *10^-15 J/s (reminder 1 watt = 1 J/s)
intensity<-1.13 *10^-15 
# energy photon = 3.94*10^-19 J
# E = hc/lambda (eV)
hc<-1.24 *10^-6 # eV/ m
hc<-1.24 # eV/micro m
E<-hc/lambda 
# one eV = 1.602 *10^-19 J 
result<-1.602 *10^-19*E  # result in J
number_photons<-intensity/result

4.2 Function for any value of wavelenght (lambda in micro m)

numberphotons<-function(lambda) {  # see PPFD function comments to understand how work the function. Input = lambda
  
  intensity<-1.13 *10^-15 
  hc<-1.24 
  E<-hc/lambda 
  denom<-1.602 *10^-19*E  
  result<-intensity/denom
  return(result)
}

4.3 Graph giving the number of photons for any value of wavelenght

lambda_values<-seq(400,800,by=50)
ggplot(data=NULL,aes(x=lambda_values,y=numberphotons(lambda_values))) + ### ggplot use a dataset (here null because I used a vector of values, aes is to define x and y)
  geom_point() + # I used geom_point because I want points in the graph, for example if you want a line you should use geom_line
  ylab("number of photons") + # specify y axis legend
  xlab("wavelength (nm)") + # specify x axis legend
  ggtitle("Dark vision: Number of photons as a function of wavelenght") # specify the title

### ggplot use a dataset (here null because I used a vector of values, aes is to define x and y)

5.Calculate CO2/O2 ratio in plant

#pick one CO2 value  (this is from Zhu's paper Fig.3)
     #present   400ppm
     CO2<-0.0004 
     #past  #220ppm,parts-per-million
     CO2<-0.00022
     #future   #700ppm
     CO2<-0.0007

#this is O2 consentration
O2<-0.21  
cat("CO2 concentration:",CO2,"%","\n")

## CO2 concentration: 7e-04 %

cat("O2 concentration:",O2,"%","\n")

## O2 concentration: 0.21 %

# step 1:get into plant (physically diffusion)
ratio1<-O2/CO2
cat("The O2/CO2 ratio of physically diffusion:",ratio1,"%","\n")

## The O2/CO2 ratio of physically diffusion: 300 %

# step 2:water dissolving
alpha<-25  #This is the CO2/O2 Solubility ratio at 25 degree
ratio2<-ratio1*1/alpha  
cat("The O2/CO2 ratio of water dissolving:",ratio2,"%","\n")

## The O2/CO2 ratio of water dissolving: 12 %

# step 3:
tau<-75  #This is the CO2/O2 Specificity ratio 
ratio3<-tau/ratio2
cat("The O2/CO2 ratio of water dissolving:",ratio3,"%","\n")

## The O2/CO2 ratio of water dissolving: 6.25 %

#results:
#past     ratio3:1.964286
#future    ratio3:6.25

6.Relationship of efficiency and temperature

# equations are from paper: Improved temperature response functions for models of Rubisco-limited photosynthesis

  
  
  ###this is the function of tau with temperature (equation 6 in the paper)
      #This function will call another function:parameter(Tk)
      #inputs: c:scaling constant (dimensionless);
      #        Ha:energy of activation (kJ mol-1);
      
      R=8.31445  #molar gas constant, M2kgs-2K-1mol-1
       Temp.f = function(c,Ha){
        parameter = function(Tk){
        #Tk is leaf temeprature
          exp(c-Ha*10^3/(R*(Tk+273.15)))
        }
        parameter
      }
  
  #This is the function of leaf temeprature
      #all c and Ha value are from table 1 in the paper
      #these four variables are functions with input Tk 
      Vc.max = Temp.f(26.35,65.33) #maximum RuBP saturated rate of carboxylation (mmol m-2 s-1);
      Vo.max = Temp.f(22.98,60.11) #maximum RuBP saturated rate of oxygenation (mmol m-2 s-1);
      Ko = Temp.f(20.3,36.38)  # Michaelis constant for O2 (mmol mol-1);
      Kc =  Temp.f(38.05,79.43 )# Michaelis constant for CO2 (mmol mol-1);

   #efficiency
        #this is from Zhu's paper
        max.conversion=0.487 #we get this value from page 154
        #input: CO2 concentration
        #       O2 concentration
        #       tempreature
        #output:Decrease value by photorespiration
        #       Efficiency

        efficiency = function (CO2,O2,temp){
          tau = Vc.max(temp)*Ko(temp)*1000/(Kc(temp)*Vo.max(temp))
          phi = O2/(CO2*tau) #this is the ratio of RuBP oxygenation to carboxylation
          dpr = 1-3*(1-0.5*phi)/(3+3.5*phi)  #decrease in  caused by photorespiration
          conversion=max.conversion-dpr#energy conversion efficiency
          list(dpr,conversion)
        }

#prepare for plot
    #temperature range
    leaf.temp=c(5:40)

    #CO2 concentration
    CO2=0.00038

    #O2 concentration
    O2=0.21

    #efficiency
    effici = efficiency(CO2,O2,leaf.temp)

    dpr=effici[[1]]
    conversion=effici[[2]]

#plot-------------------------------------

    #plot four intermediate variables

    plot(leaf.temp,Vc.max(leaf.temp),type="l",xlab="Temperature",ylab="Maximum RuBP saturated rate of carboxylation (mmol m-2 s-1)")

    plot(leaf.temp,Vo.max(leaf.temp),type="l",xlab="Temperature",ylab="Maximum RuBP saturated rate of oxygenation (mmol m-2 s-1)")

    plot(leaf.temp,Ko(leaf.temp),type="l",xlab="Temperature",ylab="Michaelis constant for O2 (mmol mol-1)")

    plot(leaf.temp,Kc(leaf.temp),type="l",xlab="Temperature",ylab="Michaelis constant for CO2 (mmol mol-1)")

    #Decrease value by photorespiration channges with temperature
    #dimensionless
    plot(leaf.temp,dpr,type="l",xlab="Temperature",ylab="Decrease by Photorespiration")

    #energy conversion efficiency channges with temperature
    #dimensionless
    plot(leaf.temp,conversion,type="l",xlab="Temperature",ylab="Energy Conversion Efficiency")

7.Device - Read data

library(readr)
#tab<-read_csv("https://raw.githubusercontent.com/Agron508/Homework/master/SS-110_1010_12918noon.csv")
tab<-read.csv("https://raw.githubusercontent.com/Agron508/Homework/master/SS-110_1010_12918noon.csv")
library(ggplot2)

Some cleaning

data<-tab[2:482,] # select rows
colnames(data)[1]<-"waveL"
data$waveL<-as.numeric(as.vector(data$waveL)) # to get a x_continous scale, data need to be numeric not factor

Plot the Irradiance as a function of wavelenght (nm)

time<-data[,84]# here I just pick a column randomly
class(time)

## [1] "factor"

# time has to be converted into numeric 
as.numeric.factor <- function(x) {as.numeric(levels(x))[x]} # function to convert all the level of the factor into numeric without losing information
time<-as.numeric.factor(time)

## Warning in as.numeric.factor(time): NAs introduced by coercion

ggplot(data=data,aes(x=waveL,y=time))+  # see part 4.3 to find some comments about ggplot
  geom_line() +
  scale_x_continuous(breaks = round(seq(min(data$waveL), max(data$waveL), by = 100),1)) + 
  scale_y_continuous(breaks = round(seq(min(time), max(time), by = 0.005),1)) +
  xlab("wavelength (nm)") +
  ylab(expression(paste("Spectral Irradiance (W/m2/",mu,"m)")))

8.Calculate the biomass (carbon)

St = total incident solar radiation during growing season of which 48.7% is photosynthetically active (MJ/m2) Ei = light interception efficiency (between 50-75%) Ec = energy conversion efficiency Ep = partitioning efficiency or harvest index (here ignored) Yp = yield potential (MJ/m2)

Seed_ec = seed energy content (MJ/kg) = 22.2 in average Leaf_ec = leaf energy content (MJ/kg) = 19 in average Stem_ec = stem energy content (MJ/kg) = 17 in average Total_ec = average energy content for the plant = (22.2+17+19)/3=19.4 (MJ/kg)

St=2000
Ei=0.6 
Ec = 0.06 # at 35°C
Yp = St*0.478*Ei*Ec

Total_ec=19.4 

biomass=Yp/Total_ec # biomass (carbon) of the whole plant 

# the biomass is equal to 1.77 kg/m2

8.1 3D plot of biomass

x<-c(0.5,0.6,0.7,0.75) # x = Ei
y<-c(0.05,0.1,0.2,0.3) # y = Ec
St<-2000
data<-expand.grid(x,y) # gives all the combinations of x and y
colnames(data)<-c("x","y")

data$z<-St*0.478*data$x*data$y  # biomass in MJ/m2

library("plot3D")

scatter3D(data$x, data$y, data$z, phi = 0, bty = "g",
          pch = 20, cex = 1.8, 
          ticktype = "detailed",
          xlab = "light interception efficiency",
          ylab ="energy conversion efficiency", 
          zlab = "yield potential (MJ/m2)",
          main="Yield potential (MJ/m2), Ei and Ec")

library(plotly)

## 
## Attaching package: 'plotly'

## The following object is masked from 'package:ggplot2':
## 
##     last_plot

## The following object is masked from 'package:stats':
## 
##     filter

## The following object is masked from 'package:graphics':
## 
##     layout

x<-c(0.5,0.6,0.7,0.75) # x = Ei
y<-c(0.05,0.1,0.2,0.3) # y = Ec
St<-2000
data<-expand.grid(x,y) # gives all the combinations of x and y
colnames(data)<-c("x","y")
data$z<-St*0.478*data$x*data$y  # biomass in MJ/m2

p <- plot_ly(data, x = ~x, y = ~y, z = ~z, color = ~z) %>%
  add_markers() %>%
  layout(scene = list(xaxis = list(title = 'Ei'),
                      yaxis = list(title = 'Ec'),
                      zaxis = list(title = 'Yp')))
p

9. Carbon assimilation model

## need to find a way to calculate C

cass <- function(result){
  
  I <- 0.035   ## with the assump. that r = 3 and O = .21 equation (3)
  vc <- 98    ## constant taken form Fig 2 explaination
  C <- 270    ## micromol per mol
  Rd <- result  
  
  ## A = (1-I/C)vc-rd
  
  A <- (1 -(I/C))*vc -Rd
  
  return(A)
}